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Abstract 

We present a model for the motion of an average atom in a liquid 
or supercooled liquid state and apply it to calculations of the velocity 
autocorrelation function Z(t) and diffusion coefficient D. The model 
trajectory consists of oscillations at a distribution of frequencies char- 
acteristic of the normal modes of a single potential valley, interspersed 
with position- and velocity-conserving transits to similar adjacent val- 
leys. The resulting predictions for Z(t) and D agree remarkably well 
with MD simulations of Na at up to almost three times its melting 
temperature. Two independent processes in the model relax velocity 
autocorrelations: (a) dephasing due to the presence of many frequency 
components, which operates at all temperatures but which produces 
no diffusion, and (b) the transit process, which increases with increas- 
ing temperature and which produces diffusion. Because the model 
provides a single-atom trajectory in real space and time, including 
transits, it may be used to calculate all single-atom correlation func- 
tions. 

1 Introduction 

In order to explain the experimental fact that the specific heat of a solid 
changes little, while its self-diffusion coefficient changes greatly, when it melts 
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to a liquid, Frenkel JT|, suggested that an atom in a liquid undergoes ap- 
proximately harmonic vibrations about an equilibrium position, occasion- 
ally jumping from one equilibrium position to another; these jumps are re- 
sponsible for self-diffusion. Using molecular dynamics calculations, Stillinger 
and Weber ([§-0) demonstrated the presence of local many-particle minima 
("inherent structures") in the potential surface underlying the liquid state, 
and they observed that "diffusion and fluid flow within a liquid may be in- 
terpreted as transitions between . . . local minima" ||. Building on these 
ideas, Zwanzig @ studied the self-diffusion coefficient D, given in terms of 
the velocity autocorrelation function Z(t) = ^(v(t)-v(0)) by the Green-Kubo 
formula || 

D= / Z(t)dt. (1) 
Jo 

For harmonic motion about a many-particle equilibrium position, Z(t) is 
given by 

kT r 

m= ~M y^) c ° s M^ ( 2 ) 

where p[iS) is the density of normal mode frequencies. (The derivation of 
this formula is discussed in |7]] and performed in Section II of ||31|| ■) Zwanzig 
suggested that jumps between equilibrium positions will have the effect of 
multiplying this expression by a factor exp(— t/r zw ), where the "hopping 
time" r zw is characteristic of the time between jumps. Much effort has been 
devoted to developing these ideas into full harmonic theories of liquid dynam- 
ics, particularly theories of self-diffusion in liquids and supercooled liquids 
dPI-PH). In most theories one finds p{ui) by expanding the potential energy 
to second order around each of some set of configurations, diagonalizing the 
second-order term in the potential (often called the dynamical matrix) to find 
the frequency distribution for that configuration, and averaging over all con- 
figurations chosen. This is done in one of two different ways. In quenched 
normal mode (QNM) theories || £10], [II], £I3|, [JOJ, p{u) is calculated at 
several potential minima and averaged, while in instantaneous normal mode 
(INM) theories (0-0, @, 0) p(u) is averaged over a thermal 

distribution of configurations, with no special emphasis placed upon configu- 
rations in potential valleys. This difference manifests itself in the fact that in 
INM theories, the configuration-averaged p(u>) usually includes both real and 
imaginary frequencies (corresponding to stable and unstable normal modes), 
since most configurations will not lie at the bottoms of valleys, while in QNM 
theories only real frequencies are represented. In addition, the QNM den- 



2 



sity of states for a given system will in general be temperature-independent, 
while the INM density will depend strongly on temperature. The two most 
prominent ways to determine t zw are (a ) to extract it from the imaginary fre- 
quency INM distribution, developed most notably by Keyes ([K|-|2^]), and 



(b) to set equal to the long-time decay rate of the "cage correlation func- 
tion" of Rabani, Gezelter, and Berne ^2|. (Another theory uses Cao and 
Voth's frequency-dependent multiplicative factor |44], [4^].) Notice that none 
of these theories attempt to model the actual motion of a diffusing particle 
in the liquid and then calculate Z(t) from this motion; they try to model the 
effects of diffusion on the autocorrelation function directly. 

The theory of Z(t) we propose differs from those discussed above both 
in how we determine p(u) and how we model the effects of diffusion. In 
Section 0, we suggest a density of states of the QNM type that also incor- 
porates insights drawn from studies of the potential energy surface of liquid 
Na by Clements and Wallace EB. These studies were undertaken to test 



a theory of monatomic liquid dynamics proposed by Wallace [52| that has 
been applied previously with some success to the thermodynamics of a wide 
variety of liquid metals |53| and an earlier study of Z{t) |54| , and the present 
work is also intended to lend credence to that theory. We propose a model 
for the process of diffusion based on the following observations. The sys- 
tem moves in a set of nearly harmonic many-particle valleys, and the motion 
within each valley may be analyzed into normal modes of vibration about 
the valley minimum. The motion of the system from one valley to another 
is called a transit, and it corresponds to a change in the equilibrium posi- 
tions of a small group of atoms. When a transit occurs involving a given 
atom, or one of its neighbors, the normal mode eigenvector components for 
that atom will change, and since this can happen many times during a single 
vibrational period (as we will see in Section ^), the eigenvectors will not nec- 
essarily provide a useful basis for describing the motion. This suggests that 



an independent atom model will be a good theoretical starting point |54 
(Notice that this argument does not apply to INM, where the eigenvectors 
change continuously with the motion, instead of discontinuously only at tran- 
sits.) Therefore, we propose a mean-atom-trajectory model which describes 
a single average particle in the diffusing liquid actually transiting between 
single-particle equilibrium positions, and from this model we calculate Z(t) 
directly. In Section |3| we compare our predictions with molecular dynamics 
(MD) simulations of liquid Na over a very broad range of temperatures, and 
in Section II we discuss our results. 
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2 The model 



2.1 Density of states 



Wallace |5lJ has predicted that in any monatomic liquid the many-body 
potential valleys can be divided into three categories: the few crystalline val- 
leys; so-called "symmetric" valleys which retain some remnant of the crystal 
system's symmetry; and "random" valleys, in which no crystal symmetries 
remain. He argued further that the random valleys should greatly outnum- 
ber the symmetric ones (thus controlling the statistical mechanics of the 
liquid), and that all random valleys should have the same distribution of 



normal mode frequencies. Wallace and Clements |55|, [56| have verified these 
predictions for liquid Na, as well as discovering criteria that one can use to 
determine whether the system is in a symmetric or random valley when non- 



diffusing. From Fig. 7 of |55| one can construct a distribution p{ui) for liquid 
Na that will be valid whenever the system remains in a random valley; this 
p{uj) is shown in Fig. [[} However, since we actually have the set of normal 
mode frequencies found in J55|] , we can use them in Eq. (^) directly, so 
our nondiffusing Z{t) is 

1 kT 

where N is the number of particles in the system, and the number of normal 
modes is 3iV — 3 because the three zero-frequency modes corresponding to 
center of mass motion are not excited. This is our model Z(t) when the 
system remains in a single random valley without diffusion. Note that be- 
cause all of the random valleys have the same p(oo), an average over quenched 
configurations is unnecessary. 



2.2 Motion of an Average Particle 

Our first goal is to construct a model for the motion of an average particle 
that reproduces Eq. (|3|) for Z{t) in the absence of transits. Since the system 
is transiting with overwhelming likelihood from random valley to random 
valley, and since all the random valleys have the same frequency distribution, 
it is sensible to model an average particle's motion in terms of oscillations at 
those frequencies, or 

r{t) = R + u(t) 
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Figure 1: p(u) for liquid Na constructed from the set of frequencies 
7 of f55|. Here 5t is the MD timestep (1.4 x 1(T 15 s). 



= R + ^2w\sm(u\t+ ax), (4) 

A 

where the particle's position r(t) is divided into a center of oscillation R and 
oscillations u(t) about that center, and the parameters in u(t) aside from the 
u\ have yet to be determined. Let us assume that the values of the phases a>\ 
are randomly distributed among the particles; then one calculates Z(t) from 
Eq. (|j) by differentiating to find v(t), computing the product v(t) -v(0), and 
averaging over each of the a\ separately; the result is 

Z ( t ) = r \ w *\ 2u l cos(u x t). (5) 
b A 

Eq. (|]) becomes Eq. (|3|) with the choice 



_L 2kT A 

where W\ is an arbitrarily chosen unit vector. Thus Eq. (Q) with the phases 
a\ randomly chosen and W\ given by Eq. (^j, with the unit vectors W\ 
also randomly chosen, constitute our mean-atom-trajectory model when the 
system is not diffusing. 

To include diffusion in our model, we must incorporate both the rate at 
which transits occur and their effect on the particle's motion. We will allow 
the rate v to be a temperature-dependent parameter that we will determine in 
Section ^| by fitting to MD simulations (so the probability of a transit in small 
time At is uAt), and we assume that the transit occurs instantaneously (the 
particle simply crosses the surface separating distinct valleys), so it must 
conserve both the particle's position r(t) and velocity v(t). To be more 
specific, we assume that the transit occurs in the forward direction, so that 
the center of the new valley lies an equal distance away from the particle 
but on the opposite side from the center of the old valley. Let r ore (t), 
i? bcfore , and u hcime (t) be the position parameters from Eq. (f|) before the 
transit, and let r after (t), i? aftcr , and n aftcr (t) be the parameters after; then 
our assumption of a forward transit implies that u r (t) = —u hetore (t), and 
this together with r before (t) = r after (t) implies 

R aiter = R heioie + 2n before (t). (7) 

This is the change in R produced by a transit. We choose to leave the 
unit vectors ib\ in Eq. (B|) unaffected by transits, leaving only the effect 
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on the phases a\ to be determined. They must change in such a way as 
to reverse the sign of u(t) but conserve v(t); since u(t) is a sum of sines 
while v(t) is a sum of cosines, this is easily done by reversing the signs of 
the arguments (u\t + a A ) in Eq. (f|). Let the transit occur at time t Q ; then 

coxt + af er = -(u x to + a h x ciorc ) so 

af er = -2u x t - a h x efore . (8) 

Thus, a transit is implemented at time t by leaving the Wx alone and making 
the substitutions 

R -> ii + 2u(t ) 

a A -> -2u; A t - a A . (9) 

This conserves r(£), reverses the sign of and conserves v(t). 

Now the model consists of two parts, (a) Between transits, the average 
particle moves nondiffusively as given by Eq. (|J) and (^j), with the phases a A 
and unit vectors wx assigned randomly, (b) In each small time interval At 
a transit occurs with probability uAt; if it occurs, it replaces R and the a A 
with new values according to Eq. (^[). With the addition of transits, we can 
no longer express r(t) and v(t) in closed form at all times, so we no longer 
have a closed form for Z(t); but the model can be implemented easily on a 
computer, and then the data from the run can be used to calculate Z(t) in 
a manner analogous to an MD simulation. We turn to comparison of the 
predictions of this model with MD results next. 



3 Comparison with MD 

The MD setup used to test our model is that described in |5l| with two 
changes: iV = 500 in all runs and the MD timestep was reduced to 8t = 0.2t*, 
where t* = 7.00 x 10~ 15 s is the natural timescale defined in [J55 ]. (The 



system's mean vibrational period r = 27r/co> rms , where u; rms is the rms average 
of the normal mode frequency distribution, is approximately 287 St.) We 
performed equilibrium runs of the system at 6.69 K, 22.3 K, 216.3 K, 309.7 
K, 425. OK, 664.7 K, and 1022.0 K; at the lower two temperatures the system 
is not diffusing (as can be seen from the system's mean square displacement, 
or from the integral of Z(t)), and the system is diffusing otherwise. Since 
T m = 371.0 K for Na at this density, our simulations range from the glassy 
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regime to nearly three times the melting temperature. We then ran the 
model for various values of v, adjusting until the model matched the value of 
the first minimum of Z(t) at each temperature. The values of v that we fit 
for all temperatures are given below; Figs. through [7] compare the model's 
predictions with the MD results for Z(t) = Z(t)/Z(0). The model requires 
v = for both nondiffusing states, so they are presented together in Fig. |^. 

T(K) v(r-i) 
6.69 0.0 
22.3 0.0 
216.3 0.35018 
309.7 0.60276 
425.0 0.83985 
664.7 1.24858 
1022.0 1.68774 

Notice that in all diffusing cases v is of the same order of magnitude as r -1 , 
indicating roughly one transit per mean vibrational period, as predicted in 
[ 57j and noted in Section |l[ 

The most obvious trend in Z(t) is that its first minimum is rising with 
increasing T; this is the primary reason for the increasing diffusion coefficient 
D. Note that the model is able to reproduce this most important feature 
quite satisfactorily. In fact, all fits of the model to the MD results capture 
their essential features, but we do see systematic trends in the discrepancies. 
First, note that the location of the first minimum barely changes at all in 
the model as v is raised, but in MD the first minimum moves steadily to 
earlier times as the temperature rises. The first minimum occurs at a time 
roughly equal to half of the mean vibrational period (recall r = 287 5t), 
so the steady drift backward suggests that the MD system is sampling a 
higher range of frequencies at higher T. Also, for the three lowest diffusing 
temperatures the model tends to overshoot the MD result in the vicinity of 
the first two maxima after the origin, and at the highest two temperatures 
this overshoot is accompanied by a positive tail that is slightly higher than 
the (still somewhat long) tail predicted by MD. These overshoots should 
clearly affect the diffusion coefficient D. To check this, we calculated the 
reduced diffusion coefficient D, the integral of Z(t), which is related to D by 
D = (kT/M)D. The results are compared to the values of D calculated from 
the MD runs in Fig. [8[ In all of the diffusing cases, the model overestimates D 
by roughly the same amount, which we take to be the effect of the overshoots 
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Figure 3: The model prediction for Z{t) at v = 0.35018 r 1 compared with 
the MD result for supercooled liquid Na at T = 216.3 K. 
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Figure 4: The model prediction for Z{t) at v = 0.60276 t~ 1 compared with 
the MD result for supercooled liquid Na at T = 309.7 K. 
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Figure 5: The model prediction for Z(t) at v 
the MD result for liquid Na at T = 425.0 K. 
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Figure 6: The model prediction for Z(t) at v 
the MD result for liquid Na at T = 664.7 K. 



I 

800 1000 



1.24858 t 1 compared with 



13 




-0.2 - 



-0.4 



200 400 600 800 1000 
t (St) 



Figure 7: The model prediction for Z(t) at v — 1.68774 r 1 compared with 
the MD result for liquid Na at T = 1022.0 K. 
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Figure 8: D as a function of T for both the model and MD. 
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at the first two maxima. At the higher temperatures the discrepancy is also 
higher, presumably due to the model's long tail. 



4 Conclusions 

We have presented a single-atom model of a monatomic liquid that provides a 
unified account of diffusing and nondiffusing behavior. The nondiffusing mo- 
tion is modeled as a sum of oscillations at the normal mode frequencies (Eq. 
(|4]) and self-diffusion is accounted for in terms of instantaneous tran- 
sits between wells, which occur at a temperature-dependent rate v. Since 
this model gives a simple and straightforward account of the motion itself, 
it can easily be used to calculate any single-atom correlation function one 
wishes; here we have focussed on Z(t) and its integral D. The relaxation 
of correlations expressed by the decay of Z(t) arises here from two distinct 
processes: Dephasing as a result of the large number of frequencies in the 
single-valley motion, and transits between valleys. The dephasing effect pro- 
duces relaxation but not diffusion: It causes Z(t) to decay but its integral 
remains zero. On the other hand, transits certainly contribute to relaxation, 
but in addition they raise the first minimum of Z(t) substantially, increasing 
its integral and providing a nonzero D. 

Most other workers in the field have studied Lennard- Jones or molecular 
liquids; the only other work with liquid Na we have found is Wu and Tsay's 



INM analysis jig, |49[], with which our results are of comparable quality. 
This is remarkable in light of the fact that our model of the transit process 
is exceedingly simple; one would expect that a more realistic model would 
produce even better results. We are also pleased to see that the model 
retains its validity from the glassy regime to well beyond the liquid's melting 
temperature. 

In light of these results, answers to the following questions are worth 
pursuing. Do other monatomic liquids exhibit the same division of their po- 
tential valleys into random and symmetric that liquid Na does? (It is known 



that LJ Ar does |58|.) How harmonic are these valleys? How similar are 
their frequency distributions? Fig. shows that in Na the valleys are very 
nearly perfectly harmonic, and we expect the same qualitative potential sur- 
face for all nearly-free-electron metals (a total of 24 elemental liquid metals), 
but other liquids might show more pronounced anharmonicities, and those 
would need to be accounted for in the model. Could a more sophisticated 
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model of the transit process (as opposed to simply transiting forward) pro- 
duce the shift in the first minimum and smaller long-time tail shown in MD? 
Can one develop a theory to predict the transit probability vl (Such a theory 
would be conceptually related to the decay of the cage correlation function 
of pOj 62|.) Finally, how can these ideas be applied to theories of other 



transport coefficients, such as bulk and shear viscosities? Future work will 
focus on answering these questions. 
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